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ABSTRACT 


A computer program based on the dual reciprocity 
boundary element method has been developed for solution of 
nonlinear transient heat conduction problems in two- 
dimensional homogeneous isotropic regions. In the present 
work, nonlinear boundary conditions have also been incor- 
porated. Further, the concept of the residual of input- 
output thermal energy balance has been used to assess the 
accuracy of computed ' solutions . This is particularly 
useful for nonlinear problems. Its application to boundary 
element solutions in the present work attests to the accuracy 
of the computed solutions. For problems solved here, the 
energy-residual was found to be of the order 0.5% or less. 

A least-squares time integration scheme was introduced, 
and was found to be computationally more efficient and 
accurate than other schemes such as the fully implicit, the 
Crank-Nicholson and the Galerkin schemes. Numerical experi- 
ments on convergence indicate second order uniform convergence 
with spatial mesh-refinement for the fully implicit scheme. 

The rate of convergence for this scheme is of first order 
with respect to time increment. With the least-squares scheme, 
the convergence of the boundary solutions with spatial mesh- 
refinement is not uniform, whereas, quadratic (or higher) order 
convergence with respect to time increment is observed. 



CHAPTER 1 


INTRODUCTION 

1 . 1 Motivation 

Nature is inherently nonlinear. Hence, the 
of 

occurrence/nonlinear processes is a rule rather than an 
exception. It is the nature of nonlinearities present 
in a process which dictates whether we can approximate 
it by a linear model or not. If the nonlinearities are 
weak, we can obtain significant insights into the process 
by analysing its linear model. Otherwise, we have no 
option except to model the process as it exists, with 
suitable approximations. Nonlinear heat conduction 

problems belong to this later class. These are encoun- 
tered in a variety of processes of engineering interest 
such as thermal and nuclear power generation, space-appli- 
cations, manufacturing processes etc. 

Nonlinearities in a heat conduction process may 
arise from: 

(i) Temperature dependent thermal diffusivity. 

(ii) Nonlinear boundary conditions, e.g. due to 
heat radiation. 

(iii) Nonlinear heat sources. 

(iv) Moving interface, e.g. due to phase change. 
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The first two types of nonlinearities arise mostly 
in high-temperature applications, and these are addressed 
in the present work. The motivation has obviously come 
from a wide range of problems already indicated above. 

1 . 2 Reasons for Interest in Boundary Element Techniques 

Mathematical model of a transient heat conduction 
process gives us an initial-boundary value problem. This 
problem is, in general, not amenable to closed form analytical 
solutions. This is true for nonlinear problems in parti- 
cular. Even for a linear problem which admits an analytical 
solution, the result is often in the form of integrals or 
infinite series. The above features have been the main 
reasons for interest in approximate numerical solution 
procedures . 

The modern numerical techniques have been, made possible 
by phenomenal advancements in digital computer-technology in 
the past three decades. The finite element and the finite 
difference methods have been the two most popular methods. 

The application of boundary element techniques to the transient 
heat -conduct ion problems has been a comparatively recent 
phenomenon . 

The interest in the boundary element techniques stems 
from the fact that it requires, for the most part, the discre- 
tization of the boundary only, as compared to the finite 
element or finite difference methods which require the 
discretization of the entire domain. Even for' the transient 
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problems which may require incorporation of initial conditions 
through domain integrals, the final system of equations is of 
the order of boundary unknowns only. 

In the boundary element method, the governing differ- 
ential equation is transformed into an equivalent integral 
equation using either an indirect formulation (i.e. by repre- 
senting the unknown function as a single or double layer 
potential) or a direct formulation employing a weighted 
residual statement in conjunction with free space Green's 
functions. Owing to its conceptual simplicity and wide range 
of applicability the latter approach has been the most popular. 
The integral equation so obtained is solved in the same manner 
as in the finite element methods. 

Various formulations proposed for the solution of 
transient heat conduction problems with the boundary element 
techniques include the use of Laplace transforms [l], finite- 
difference approximation of the time derivative [2], the time 
dependent fundamental solutions [3, 4 ], and the dual reciprocity 
boundary element method [5-8] . The latter offers distinct 
advantages over the others as it involves only boundary 
integrals retaining the 'boundary -only ' character of the 
method. It utilizes fundamental solution to Laplace' equation 
which is only space dependent. This formulation has been 
extended to deal with nonlinear diffusion problems with 
linear boundary conditions [8 ] . We use this formulation 
in the present work. 


Numbers in brackets denote References. 
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1 ^ Qbj ectives o-f the' Present Work 

Previous applications of the dual reciprocity 
boundary element method for linear as well as nonlinear 
conduction problems [5-8 3 have considered only linear 
boundary conditions. Time-integrations, have been performed 
using only point-collocation schemes. Also, no systematic 
study has been reported concerning the convergence of 
numerical solutions with respect to mesh-refinement and 
decreasing time-increment except a tentative study presented 
in reference [8] . We take the above tasks in the present 
work. In order to achieve any of the above, our first task 
is to develop a comprehensive computer program based on 
the chosen formulation. We restrict ourselves to two- 
dimensional problems only. The program should be capable of 
dealing with linear as well as- nonlinear conduction problems 
involving time-dependent or time-independent boundary 
conditions. We attempt to introduce the least-squares scheme 
in time for linear conduction problems and assess its 
viability in the boundary element techniques. We attempt a 
numerical study of convergence of the boundary solutions for 
the point-collocation, and the proposed least-squares time- 
integration scheme. As the final step, we propose to solve 
problems involving nonlinear boundary conditions. To assess 
the accuracy of numerical solutions, we propose to introduce 
the unbalance of input-output thermal energy as a measure of 
accuracy and physical acceptability of the computed solutions 
Specific objectives of this work can be listed t as: 
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i) Development of a comprehensive computer program 
for two dimensional transient conduction problems 
capable of dealing with material as well as 
boundary nonlinearities, time-dependent boundary 
conditions, and heat generation. 

ii) Incorporation of the least-squares time integration 
scheme. 

iii) Study of the convergence of the numerical solutions 
to the exact one with mesh-refinement and decreasing 
time-increment . 

iv) Comparison of the least-*- squares scheme with point 
collocation schemes. 

v) Some applications to problems involving nonlinear 
boundary conditions. Use of the energy residual 

to assess the accuracy of computed solutions. 

*■ * 



CHAPTER 2 


REVIEW OF LITERATURE 

2 . 1 General Background 

The modern theory of boundary integral equations finds ^ 
its origin in Fredholm's work in 1903. The first numerical 
application was made in 1917 by Trefftz. Most of the early 
work with boundary integral equations has been concerned with 
applications in potential problems. However, actual numerical 
implementations for problems of engineering interest appear 
only in 1960s with availability of fast digital computers. 

Various applications in potential theory, elastostatics and 
electrostatics appeared in this decade. Detailed account of 
these developments with related references can be found in 
the references [9-11 ]. The first application of boundary 
element techniques to transient heat conduction problem appeared 
in 197 0 [ 1 ] We review some representative works related 
to the transient heat conduction in the following sections. 

2 . 2 Boundary Element Methods' in Trans len t ._ Heat Conductio n 

In 1970, Rizzo and Shippy [l] proposed a boundary element 
formulation in conjunction with Laplace Transforms. Appli- 
cation of Laplace transform removes the time-dependence 
of the problem. The elliptic partial differential equation 
so obtained is used for boundary element formulation. The 
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problem is solved in Laplace transform space for sequence 
of real positive transform parameter, and numerical inversion 
is performed to compute the physical variables in the real 
space. Details of the method are available in references 
[1,10,12 ]. This method has found very limited application 
owing to its complexity. 

Another approach based on time-dependent fundamental 
solutions in the context of the direct method was proposed 
by Chang et al. [3] and Shaw [13]. This approach has been 
the most popular one and has been extended to deal with 
axisymmetric and 3-D problems of practical interest, which 
we review later in this section. 

A coupled boundary element-finite difference formul- 
ation was proposed by Brebbia and Walker [ 12 ] and Curran 
et al. [ 2 ] in 1980. This involves finite-difference appro- 
ximation of time derivative obtaining an elliptic equation. 
This is solved at each time step using the boundary element 
method. Owing to finite-difference approximation involved, 
this method imposes severe restrictions on time-increment 
to obtain accurate results. 

All the above mentioned approaches usually involve 
incorporation of the initial conditions through a domain 
integral which disturbs the ' boundary-only ' character of 
the technique. In the context of the method employing 
time-dependent fundamental solution, a time-stepping algorithm 
in which initial conditions are accounted through the 
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boundary integrals was proposed by Wrobel and Brebbia [ 1 4 3 - 
This technique, though mathematically elegant, becomes 
computationally inefficient if the number of time-steps are 
large. 

The latest approach in this series is the dual 
reciprocity boundary element method. This approach was 
proposed by Brebbia and Nardini [l5 ] in context of elasto- 
dynamics, and has been applied to the transient conduction 
problems by Wrobel et al. [5-8] . This involves approxi- 
mation of the temporal derivative by a specific type of 
product-approximation in conjunction with the fundamental 
solution of the Laplace's equation. By exploiting the 
special nature of space dependent functions chosen in above 
approximation in conjunction with the reciprocity principles, 
the domain integral is transformed into boundary integrals. 

We review this method in detail in next chapter. We refer 
to this method as DRBEM, and to the one employing time- 
dependent fundamental solutions as BEM in the following 
lines. 

Wrobel and Brebbia [l6 ] presented extension of BEM 
to deal with axisymmetric problems. They employed integral 
in ©-direction, of the three dimensional fundamental solution, 
over generating area and boundary contour of the axisymmetric 
body. This was referred to as axisymmetric fundamental 
solution. Reference [lO ] presents a detailed account of 


this method. 
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The BEM was applied to melting and solidification 
problems by Banerjee and Shaw [ 17 ] . They developed general 
boundary integral formulations for general one-, two- 
and three-dimensional problems of melting and solidification 
of metals as well as freezing and thawing of saturated 
soils, and presented their numerical solutions. 

Problems involving phase change and solidification 
were also discussed by Hong et al. [18], and Wrobel [ 19 1 

Pina and Fernandes [20 ] presented solution of 
three-dimensional transient conduction problems with harmonic 
initial conditions. The assumption of harmonic initial 
conditions allows the conversion of corresponding domain 
integral into boundary integrals, thus, requiring discre- 
tization of boundary only. Examples of heat conduction in 
a cube and unit sphere with the Dirichlet boundary conditions 
and zero initial temperatures have been presented. 

Kihara et al. [ 21 ] discussed the accuracy of boundary 

element solution for two-dimensional problems using the 

2 

stability parameters a At/ (Ax) for constant boundary 
element discretization. They also discussed the effect of 
changes in mesh pattern, based on heat flux density, on the 
numerical solution. 

Uniform convergence of boundary element solution of 
heat equations have been investigated by Onishi and Kuroki 
for pure Neumann conditions. Iso [ 22 ] extended the above 
investigation for problems with linear radiation boundary 
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conditions. Iso et al. [23] present mathematical theorems 
on uniform convergence of boundary solution and stability 
of the computing scheme for two-dimensional isotropic heat 
conduction problem involving non-isothermal boundary 
conditions. They also present numerical investigations of 
convergence of boundary solutions for mixed type problems. 
Problems involving singularities such as re-entrant corners, 
slit boundary and interzonal singularity, and infinite 
domains are also discussed. 

2 . 3 Application in Nonlinear Problems involving Material 
Nonlinearities and Nonlinear Boundary Conditions 

Koizumi et al. [24] discussed the application of 
the boundary element method to three-dimensional heat condu- 
ction problems with nonlinear boundary conditions. Onishi 
and Kuroki [25] discussed the nonlinear problems of heat 
conduction subject to radiation boundary conditions, and 
the unsteady problems subject to forced and natural 
convection using stream function and vortic'ity formulation. 

For problems involving material nonlinearities, 
Brebbia and Skerget [26 1 introduced Kirchoff’s transform 
as the new dependent variable in conjunction with constant 
diffusivity approximation and closed form variation of 
thermal conductivity. The resulting system was solved in 
Kirchoffs transform space, and analytical inversion was 
done to obtain temperature in real space. 
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Kikuta et al. [27,28] assumed linear variation of 
thermal conductivity. Using Kirchoff's transform they 
obtain a pseudo-linear integral equation in transform 
space with nonlinearities being shifted to a domain 
integral term having time derivative of transform variable 
in its kernel. The use of constant element in time, in 
conjunction with usual BEM procedure, made the evaluation 
of domain integrals corresponding to nonlinearities unne- 
cessary. Several examples were presented. 

Finally, in the context of the dual reciprocity 
boundary element formulation, Wrobel and Brebbia [8] extended 
this method to deal with nonlinear material problems by 
introducing the integral of conductivity (Kirchoff's Transform) 
as a new variable, together with a modified time variable. 

They obtained a pseudo-linear equation in transform space 
which was solved using a Newton-Raphson iteration procedure. 
Examples of problems with linear boundary conditions were 
presented. We start our work with this formulation which 
is reviewed in detail in the next chapter. 



CHAPTER 3 


DUAL RECIPROCITY BOUNDARY EL EMENT F ORMU LAT ION 

3 . 1 Mathematical Mod el of Heat Conduction Probl em 

Let q be a homogeneous isotropic heat conducting 

domain with a piecewise smooth boundary r . Let n be 

the external normal to boundary r, and, x denote the 

spatial coordinates of a field point. The transient 

heat conduction in finite time interval t < t < t„ is 

o — p 

governed by the well known equation 

Pc U = V . (KVT) + Q in a x (t t_ J . (3.1) 

dt <5 Or 

Here, P is the mass-density, c, the specific heat, K 
the thermal conductivity, T(x, t) the temperature at 
field point x at time t, Q the rate of heat generation 

g 

(per unit volume), and V denotes the gradient operator. 

The initial condition for Eq. (3.1) is given 

by 

T (x, t Q ) = T q ■ (x) on n (3.2) 

where T q (x) is a known function. 

The boundary conditions for Eq. (3.1) can, in 
general, be written as: 
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T (x, t) = T (x, t ) on r T x (t Q , ] (3.3) (a) 

q(x,t) = I (T, x, t ) on r x (t ] (3.3) (b) 

c <5 o ' r 

with r T U r q = r. In the above expressions T and 
are given functions. 

3 . 2 Non-dimensional izat ion of Governing Equations 

We define following dimensionless quantities: 


Temperature 

: u = 

(T-T ) / (T -T ) 

cL 


Spatial coordinates 

* 

: x = 

x/L 


Time (the Fourier 
Number) 

* 

: t = 

(K r /p r C r ) * {t/Ij2) 


Thermal conductivity 

★ 

: k = 

K/K r 

(3.4) 

Mass density 

* 

: p 

p/p r 


Specific Heat 

* 

: c = 

c/c r 


Volumetric heat 
generation 

q = 
y 

(Q g L 2 )/(K r (T b - 

V> 


In the above definitions, T and T. are suitably chosen 

a o 

reference temperatures. Subscript r with property values 
K, p and c indicates value of these at the chosen reference 
temperature T . L is the characteristic dimension of the 
problem. 

Introduction of definitions (3.4) in governing equations 
(3.1) - (3.3) gives us following set of non-dimensional equations 

in which we have dropped superscript 1 * ' for sake of conve- 


nience: 
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-Field Equation: 

Pc ^ = v . (kv u ) + <3 in n x (t^tp ] (3.5) 

-Initial Condition: 

u(x,0) = u q (x) on n (3.6) 


-Boundary Conditions: 


u (x, t) 

u(x,t). on r x (t ,t„] 
u o F 

(3.7) 

q (x, t) = 

f c (x,t,u) on r g X (t Q ,t F ] 

(3.8) 

with r ' u r 
u q 

= r . 



We may note here that for linear heat conduction 
problems, parameters, P , c and k in Eq. (3.5) become unity 
as per definitions (3.4). 

3 . 3 Dual Reciprocity Boundary Element Formulation for 
Linear Heat Conduction [ 5- 8-1 

3.3.1 Boundary Integral Equations [5-8 1 

From section 3.2, we have governing field equation 
for linear heat conduction given by: 

l^r = v 2 U + q g in n x (t Q , t p ] (3.9) 

As the first step in formulation of integral 
equations, we construct a relationship based on so^alled 
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reciprocity principles using the fundamental solution to 
the Laplace's equation, i.e. the generalized function u 
satisfying the equation 

- u ( 5 , x) = 6(5, x) , (3.10) 

where 6(5, x) is the Dirac delta function, and Green's 
second identity 

/ (v v 2 w - w v 2 v) dq = / (v — - w dr, 

q r 3 

(3.11) 

where V 2 is the Laplacian operator. We have from Eq.(3.11) 


* 

fin 

q 

V 2 u - u V 

2 * 
u ) dq 

* 

= fin 
r 

* 

q - q u) 

dr 

(3.12) 

with 

= ^ and 

■ 3n 

* 

q (5 

, x) 

* 

hu (f, ... 

xL 


q 

an 



Substitution 

of Eq. 

(3.10) 

into Eq,. 

(3.12) 

gives: 

/ U 

q 

where 

* 2 

v u dq = 

-c (?) 

u (?) + 

★ 

/ (u q - 
r 

* 

- q u)dr 

(3.13) 


c ( 5 ) = 

/ 6(5, 

x) dq 



(3.14) 


q 


Eq. (3.13) is the relation we shall be using in the formu- 
lation. We may note here from Eq. (3.14) and well known 
properties of the Dirac-delta distribution that if the 
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source point lies in the interior of ft , c (5 ) = 1. 

Further, if g is on the boundary r of fi , then 0 ( 5 ) 

equals the fraction of internal angle subtended by 

the boundary r at 5 , relative to the solid angle of 
ci 

the sphere in R , where d denotes the dimension of 
the problem [29] . For smooth points on r , c ( 5 ) = 1/2. 

We now construct the following weighted-residual 
statement for approximate solution of Eq. (3.9): 

/ ■ q g - V 2 u) u dfi =0 (3.15) 

Use of Eq. (3.13) in Eq. (3.15) yields: 

/ u u* dn - / q u* dn = -c(5) u(£) + / (u q-q u) dr . 

a n 9 r 

(3.16) 

in which dot stands for temporal derivative. 

For evaluation of the first domain integral in 
Eq. (3.16), we approximate u at any point in the domain fi 
by 

N . .j 

u (x, t) = If-* (x) a (t) (3.17) 

j = l 

where f-^ (x) are chosen such that there exist functions 4 ^ (x) 
satisfying the relation, 

yj _ f J _ (3.18) 
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In the above N refers to the number of functions f-^ chosen. 
With approximation (3.17), we have 

f uu an = I d J / f J u dfi (3.19) 

o j = l n 

In the light of condition given by Eq. (3.18), Eq. (3.19) 
becomes : 

• ^ * 3 o -| it 

f u u dO = z a /V 4' J u dO (3.20) 

0 j = l 0 

We now appeal to the transformation expressed by Eq. (3.13) 
obtaining 

/ uu dfl = I [ -c (E , ) ¥ J (5 ) + / (u n J 

0 j = i r 

with 

j _ 

n ~ 3n * 

Substituting Eq. (3.21) into Eq. (3.16) we finally arrive 
at the expression 

u*rt'* )dT ] d*’ 

(3.22) 

where subscript i refers to the source point g. We may note 
that the above integral equation involves only boundary 
integrals in the absence of the heat generation term. Since 


ic ic N * 

c . u • + / (q u-u q) d T = z [c + / (q - 

11 r j=i 1 1 r 


+ /q u do 


- q ) dr] ex 

(3.21) 
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the above expression has been obtained by a double application 
of reciprocity principles through transformation (3.13), 
the technique based on it is called the dual reciprocity 
boundary element method [ 5 ]. 

3.3.2 Spatial Discretization 

To model the geometry of the domain surface, we 
discretize the boundary into a series of segments r e 
which are either straight or curved lines for a 2-D problem. 
These segments are called boundary elements. Geometry of 
a boundary element is defined by the coordinates of a number 
of points inside it (called geometric nodes) and associated 
shape functions. The variations of field variables and 
other functions in a boundary element are modelled using 
interpolation functions associated with chosen points (called 
freedom nodes) in it. In this work, we have restricted our 
choice to straight line elements with constant or linear 
interpolation functions. 

While using continuous elements, the freedom nodes 
are located at the element edge. Modelling difficulties arise 
when the element edge coincides with a singularity in 
problem definition such as an abrupt change in boundary 
condition or a corner etc. At a corner, normal to the 
boundary is not defined. At an interface where the type of 
boundary conditions changes, at the same freedom node two 
different conditions are known corresponding to respective 



19 


boundary elements on the two sides. To overcome these modell- 
ing difficulties, two approaches have been used in practice. 
The first of these employs two freedom nodes close together 
on either side of problem singularity ignoring the corner 

region. The second consists in locating the freedom nodes 
in the interior of the element instead of at , the element , 

edge. The fir st • approach with double nodes at the corners, 
introduces strong linear dependence in the algebraic system 

* -j 

as functions u (g , x) and f J (x) are virtually the same at 
both nodes. Thus, - this approach is inapplicable in the 
context of the dual reciprocity boundary element method, and 
we opt for the second approach. Partially discontinuous 
elements with freedom nodes at the usual locations where 
continuity is required, and, moved into the element where a 
discontinuous interpolation is needed, have been found to be 
the best choice [ 30 ] , and we opt for these. Figure 3.1 
illustrates constant, continuous linear, and two variants 
of partially discontinuous linear elements. 

Furthermore, we divide the domain Q into N c finite 

element cells { n } . . T , with Q as a linear triangle. 

m m = 1 , N m 

c 

This division is done in such a way that the nodal points 
of the cells adjacent to the boundary coincide with the 
geometric nodes of boundary elements taken on r . Figure 
3.2 shows a sketch of discretization of the boundary and 
the domain using linear elements. 
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o Geometric nodes 
x Freedom nodes 

Fig. 3.1 Straight line elements 

(a) Constant element 

( b) Continuous linear element 

(c) Partially discontinuous linear elements. 



o 


Boundary nodes 
Internal nodes 
r e Boundary element 
S t m Finie element cell 


Fig. 3.2 Spatial discretization. 
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3.3.3 Interpolation of' Functions 

We can approximate the variation of functions u,q,y , 
and n within each boundary element using a unique set of 
interpolation functions, denoted by a column vector 0 [8]. 

We may note here that it is not necessary to approximate ¥ 
and n using interpolation .functions as these are known 
functions. However, this type of approximation brings 
significant reduction is necessary boundary integrations 
with some sacrifice in accuracy [ 8 ] . Thus, 

u| = 0 T V | r = 0 T 4^ 

e e (3.23) 

<3|r = 0 T <V i| r = f 4 

e e 

where u , q , and n J denote (column) vectors of nodal 
0 0 0 0 

values of respective functions. 

For constant elements: 0=1. The freedom node 
is usually taken at the midpoint of the element [Fig. 3 . 1 (a) 1. 
For linear elements, with two freedom nodes: 

0 = [ 0j_ 0 2 ] T (3 - 24) 

We have, for continuous linear elements [Fig.3.1(b)] , the 
interpolation functions given by [l0,30 ] 

0 X = | (1- O' 0 2 = ■ 2 U +C 1 ‘(-3.25) 
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Interpolation functions for partially discontinuous elements 
[Fig. 3. 1(c)] are [30]: 

Type 1 (Discontinuity at t= +1 Type 2_ (Discontinuity at E=- 

0 1 = 2/3 (1/2 - C) 0 1 = 2/3 Cl -0 

0 2 = 2/3 (1 +C) ' 0 2 = 2/3 (1/2 + c ) 

(3.26) 

Furthermore, we interpolate function on n using 

conventional finite element approximation. Let x k e the 

set of finite element linear interpolation functions on the 

triangle Q , then 
m 


‘g'n 


m 



(3.27) 


with 


X 



*2 



and 






T 

m 


Expressions for interpolation function x in terms of natural 
coordinates for a triangular element [Fig. 3. 3] are: 


X 


1 



x 2 " C 2 ; X 3 



(3.28) 


where 

K 3 = 1 - Cc x * 
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3.3.4 Selection of Functions ^ and n£ 

A number of possibilities, such as polynomials, 
trigonometric series, and concentrated nodal functions 
related to fixed points in space, have been suggested by 
Nardini and Brebbia [3l] . Wrobel et al. 8 used the 
set of concentrated nodal functions, referring to the 
fixed points in space as poles. We employ this set of 
functions, and in following lines we present expressions 
for these functions from reference [8 ] . 

For two- and three-dimensional problems, 


f 3 (x) = r (x, x ) 


(3.29) 


where r denotes the Euclidean distance between two points 
in space, x a field point, and Xj denotes the position of 
j th pole. Expressions for the corresponding ¥ and n are: 


, 3 (x) = ~ r 3 (x, x . ) 


n 3 (x) = ^ r 2 (x, Xj ) d (x, Xj ) 


(3.30) 

(3.31) 


with 


9 (or 12) for 2-D (or 3-D) problems, and 


d (x, x . ) = (x - x.) . n (x) 

1 


(3.32) 


n being the unit outward normal . 

Expressions for and n corresponding to a constant 

function 


f 3 (x) = . c 


(.3.33) 
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associated to an internal pole, chosen to better simulate 
the heating of the whole body by a constant value, are: 

(x) = r- r 2 (x, x.) . c (3.34) 

® J 

n"* (x) = — r ( x , x . ) d(x,x.) . c (3.35) 

b J J 

with b =4 for 2-D, and b = 6 for 3-D problems. The 
constant c is taken as maximum distance between any pair 
of boundary nodes. 

It may be noted that in order to satisfy the requirement 
of linear independence of coordinate functions f- 1 (x) , 
coincident poles (as will the case if modelling with continuous 
elements with double nodes at corners) must be avoided. 


3.3.5 The Discretized Boundary Element' Equations 

Introduction of boundary element approximations (3.23) 
in conjunction with the spatial discretization and domain 
approximations (3.27) into integral equation (3.22). yields: 


± u ± + E 6 [ ({ 0 T q* dD u e - if 0 T u*dr)q e 3 

e— 1 e ^ 


T * 


N N e 

I [c. yi. + I { (/ 0 T q*dr ) mi - a 0 T u*dr 
1 = 1 1 1 e =l r e e e 


') rig}^ a 


N T * .m 

+ I ( / X u dfi) q 

m=l U m y 


(3.36) 
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where N, N and N denote the number of poles, the number 
G c 

of boundary elements, and the number of internal finite 
element cells respectively. 

To obtain a discrete system of boundary element 
equations for nodal boundary unknowns, we employ the collo- 
cation procedure. By establishing the fulfilment of 
Eq. (3.36) at the N poles (which include all boundary freedom 
nodes plus atleast one internal pole) , and utilizing boundary 
element and internal cell connectivities, we obtain following 
system of equations: 


H U - GQ = (H¥ - Gn ) a + B Q g (3.37) 


in which U, Q, a and Q represent nodal vectors corresponding 

or 

m 

to functions u, q, a and q respectively. The elements of 

3 

matrices H, G, ¥ , n and B are constituted from the terms of 
the form: 


H. . 
U 




with 

£ = 
ir 

/ 

r j 

0 k q 

* 

(x i , x) 

dr (x) 


k 

G ir 

f , 

• 3 

0 k U 

(x i , x) 

dr (x) 


V 

ij 

V 

( x . ) ; 

l 

n ij 

= n j 


B k = 

is 

7 

* 

*k u 

(x.,x) 

do (x) 


(3.38) 
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where r = Cb(j,k) and s * Cc (m, k) . Cb and Cc denote 
global boundary element connectivity matrix and finite 
element connectivity matrix respectively. 

Evaluation of expression (3.17) at all poles obtains: 

U=Fa (3.39) 

which is used in the form 

a = F -1 U (3.40) 

Substituting Eq. (3.40) into Eq. (3.37), we get 

CU + HU - GQ - BQ =0 (3.41) 

g 

with 

C = (G n - HY) F -1 . (3.42) 

System of equations (3.41) is similar in form to the 
one obtained using the finite element method with the differ- 
ence that the presence of vector Q renders it into a ’mixed- 
type 1 system as compared to 'displacement-only 1 type finite- 
element system of equations. We can utilize any time inte- 
gration scheme for solution of this system. 
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3 .3.6 Time-Integration Schemes 

For solving the system of equations (3.41), we shall 
use a finite element discretization applied to time-domain. 

In view of infinite extent of time domain, we concentrate 
on finite domains of time, and repeat the computation for 
subsequent domains with new initial conditions [32] . In 

following lines, we assume that full domain of investigation 
corresponds with that of one element and, derive recurrence 
relations between two successive values - two point recurrence 
relations - using the weighted residual and the least squares 
formulations . 

Following standard finite element approximation proce- 
dure, we can write 

A (t ) = z (t) A j (3.43) 

j 

where A J denotes vector A at time t . , and N J (t) are scalar 
shape functions. 

Compatibility and completeness conditions require the 
N J to be at least of first order, as only the first order 
derivatives are involved in Eq. (3.41). For a typical linear 
time element shown in Fig. (3.3), interpolation and shape 
functions, and their first derivatives in terms of local 
variable 

t = . (t - t m )/At , 0 £ t <. 1 


1 3.44) 
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are : 


m 


N 


N 


m+1 


. m 

1 - z ' N 

.m+1 

C , . N 


1 

At 


(3.45 ) 


At 


with At as length of the element. 

Using the same set of interpolation functions for all 
functions viz., U, Q, and Q , in Eq. (3.41) we get: 

g 


CCH m u ra + n""* u m+1 ) + H»V + N r " +1 U ,n+1 )-G(N m Q m +N m+1 Q m+1 ) 


.m+1 


B (N m Q™ + N m+1 Q*^ +1 ) = 0 (3.46) 

g g 


A. Weighted Residual Formulation 32 

Under the assumption of one element representing the 
full domain of investigation in time, we can write a weighted 
residual statement, corresponding to Eq. (3.46), given by: 


/ W 


. rc(N m U m + n " +1 u™ +1 > + HW'V + H ra+ V +1 ) 


-G(N m Q m + N m+1 +Q ' , ’ +1 )-B(N m i nJ + N m+1 a m+1 ))dC = 0 (j = l) 

g ■ g 


(3.47) 


where represents a weighting function. 

Inserting Eq. (3.45) into Eq. (3.47), we obtain 
a recurrence relation: 


(C/At + QH)U m+1 - 9G Q m+1 = [ C/At- (1-9) H]u m + (1-9) G Q - 

+B [ (1-9) Qg +9Q™ +1 ] 


(3.48) 
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where 


1 1 
9 = / W. ? d? / / W. d? 

o o J 


Values of 9 for a series of weighting functions are: 


w . = 

J 

6(c) 


<D 

II 

O 

w, 

J 

6( C- |) 


9 = 1/2 

s; 

ii 

6 ( 1 - z) 


1 — i 

II 

a> 

w j = 

1 


9 = 1/2 

w. = 

J 

1 - c 


9 = 1/3 

w. 

J 

z 


9 = 2/3 


From stability considerations, 9 >_ 1/2 for unconditionally 
stable scheme. The fully implicit scheme corresponding to 
9 = 1 is mostly used in present work as it gives stable and 
oscillation free solution. 


B. Least Squares F or m u 1 at i o n 

In context of finite element methods in space and 
time, Zienkiewicz [32 ] has given a least-square scheme in 
time-domain, and shown that this scheme is the most accurate 
amongst various two-point integration scheme, though at a 
greater associated cost. In the following lines, we construct 
a similar scheme. Motivation for ensuing formulation is 


expected high accuracy. 
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In the context of one element representing the whole 
domain and linear interpolation functions, we write a 
functional n given by integral of square of error over the 
element, i.e. 

n = / [ c (rAj™ + N m+1 u rn+1 ) + h (N Tn u m +N m+1 u m+1 ) 

O 

- G(nV + S Xa+1 Q ,n+1 ) - B(lV + N m+1 6 m+1 )] T 

g g 

[ C (N m U m + ] dC (3.50) 

m+ 

We minimize functional n with respect to variables U m , and 
Q m+ '*' subject to prescribed constraints in form of boundary 
conditions. We can adopt two approaches. The first consists 
in introduction of a set of Lagrange multipliers, whereas 
the second involves a direct substitution process eliminating 
some of variables (equal to number of constraints) using 
constraint conditions. The first approach, though very general, 
would entail the solution of an algebraic system of order 
3N as compared to the one obtained using second approach having 
order N, N being the number of freedom nodes. Simple form 
of boundary (constraint) conditions offers further incentive 
for choice of second approach, and we opt for it. 

Introducing boundary conditions (3.7) and (3.8) and 
Eq. (3.45) into Eq. (3.50), and minimizing the functional 
with respect to vector of unknowns X defined as 
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X . 
J 



= u m+1 
J 


if boundary condition (3.7) specified at 
node j . 

if boundary condition (3.8) specified at 
node j or it is an internal point. 


We obtain following set of recurrence relations: 


— If X i = U? , then, i n equation has the form: 


_T_, 
E [{ S_C 


(c t h + h t c ) h t h 


T 

r 


(At) 


2 At 


+ 


T , G"C + G'H) yTTi+1 

'"' 1 J -t A n 


A ( 


2 A t 3 i j j 


T m+1 C^G H^G T G^G m 

{ 2At + 3 " A 3 } ij Q j " <2 At 6 A 6 >ij y j 


+ { 


T 

C C 


( At) 


(C T H - H T C) H^H _ T ( G^C _ IV n 

2 At 6 A 1 At 6 /j ij J 


m 


C T . 

11 


2 At 


u^ 1 q m T T 

(s m + s ra+1 ) + -iJ- fc 1 + s'"* 1 ) - n 

J j 5 /. j - 0 


s m 


+ S™ +1 ) } ] = o 


(3.51) C a ) 


- If X. 


q^ + 1 , then, i th equation has the form: 


_T_ 
/G C 


E tV 1 At + 3 ‘ J i j j 


n^ + qT +1 

* J T n i O J 


- { 


T 

G H 


j. . u m + 

2 At l j J 


tgTG) li -m 

c J 


+ 


rn 

(G hi 8 

n '■ 


m 


+ S 


T +1 ) 1 


(3.51) (b) 


3 


2 


0 
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In the above equation, A is a diagonal matrix, and S is a 
column vector. These are defined as: 







(Q g>k 


(3.52) 


Equations (3.51) represent the recurrence relation 
obtained from least squares formulation. Obviously this 
scheme involves more computations as compared to the scheme 
given by Eq. (3.48). 

System of equations (3.48) or (3.51) can be solved 
by using either a direct elimination procedure such as Gauss 
elimination, if the system is linear, or, a Newton-Raphson 
procedure, if it is nonlinear. 

3.3.7 Solution at Internal Points 

After we have obtained boundary unknowns and temper- 
ature^) at internal pole(s), we appeal to Eq. (3.36) to 
compute value of temperature at any specified internal 
point . Noting that for an internal point i, c_^ = 
can rewrite Eq. (3.36) as: 


1, we 
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N 

0 

u. = Z i(f r 0 T u*dr) q - (/ 0 T q*dr ) u e ] 
1 e=l *e e r e 


+ 


l [f x T u* df2) q™ 
m=l n m g 


+ 


£ [ + £ l (/_ 0 T q*dr) <p^ - (/_ 0 T u*dr) ni}]a J ’ 

j = l 1 e=l 1 e e 1 e 


(3.53) 


which after application of element connectivities takes the 
form: 


N, 


N 


N. 


b 


u i = E 

X 


E (G..Q. - H..U.) + E [¥..+ . (H. v -G.. ru .)] 

i=1 lj J ij J j = 1 L xj k=l lk kj ik kj' J 


n t 

+ Jl ^ 


(3.54) 


where N^, N and N T denote the number of boundary nodes, the 
number of poles, and the number of total nodes (boundary + 
internal) respectively. Values of can be calculated using 

Eq. (3.40), and interpolation given by Eq. (3.45), i.e., , 


U, 


N -1 

■ -t < 1 ~ < 


(3.55) 


and 
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Fluxes at internal points can be computed by taking derivative 
of Eq. (3.53) with respect to coordinates of the internal 
point . 

3 . 4 DRBEM Formulation for Nonlinear Conduction f'8 ] 

The general heat conduction equation can be written 

as: 


pc = v. (k V u) + q (3.5) 

o'- • y 

with associated initial and boundary conditions given by 
Eqs. (3.6) - (3.8). In order to obtain a boundary element 

formulation corresponding to the above problem, Wrobel 
and Brebbia [ 8 ] introduced the integral of conductivity as 
new dependent variable in conjunction with a modified time 
variable obtaining a linear conduction equation in Kirchhoff 
transform space. In following lines, we review the above 
mentioned formulation and the solution procedure based on 
it. 

3.4.1 Transformation of' Governing' Equations 

We introduce a new dependent variable defined as: 

u 

U = / k (u) du (. 3.56) 

u r 

where u denotes an arbitrary reference value. The above 
r 
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relation represents the so-called Kirchhof f ' s transformation. 
We denote this transformation rule by T, i.e. 


U = T (u) 


(3.57) 


From equation (3.56), we have 


VU = k v u 


(3.68) 


3 U , 3 u 

at K at 


(3.59) 


Introducing Eqs. (3.58) and (3.59) into (3.5) obtains: 


1 81J 

a at 


v 2 u 


(3 .60) 


with a = a (u) = k/pc, the non-dimensional thermal 

diffusivity. Since u is a continuous function of x and t, 
we can write a = a (x,t) . We now define a modified time 
variable x by the relation: 


t 

x = f a (x, t) dt 
t 

o 

From Eq. (3.61), we get 


3 x 
at 


a 


(3.61) 


(3.62) 


Substituting Eq. (3.62) into Eq. (3.60), we obtain 
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3 U 2 * 

i7 = v U + q g in fi x (t Q , t p ] (3.63) 

Initial and boundary conditions in the Kirchhoff Transform 
space are given by: 

U(x, 0) = T (u (x) ) on ft (3.64) 

U (x, t) 

Q (x, .t) 

We note that the Eq. (3.63) has the same form as. the 
linear conduction equation (3.9) with the exception that the 
modified time variable in it is itself a function of x. 

This necessiates the employment of an iterative solution 
process, which is outlined in next sub-section. 

3.4.2 Solution Process 

The application of the dual reciprocity boundary 
element method to Eq. (3.63) obtains a system of equations 
similar to Eq. (3.48) with the weighted residual formula- 
tion. The recurrence relation can be written as: 


= U (x, t ) = T (u (x, t ) ) on r x (t ,t„] 

U Or 

= f~ = q(x,t) = f c (x, t, T~ 1 (U) ) 

on r x (t , t„l 
q o FJ 
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(£ + ©H) U m+1 - ©G Q m+1 =[c - (1-9) H ] U m + (1-0)G Q m 

+ B[ (1—0) + © Qg +1 ] (3.66) 

with C. . = C. ./at. . In above equation U and Q 

1J iJ J 

represent nodal vectors in the Kirchhoff Transform space, 
and Atj denotes the step value of the modified time 
variable t at node j, i.e.. 

At • = a , At (3.67) 

J J 

For the special case of constant diffusivity, i.e., 
a = 1, Ax j = At for all points j. In this case, we 
can employ a standard equation solver to obtain the 
solution in transform space. For general case, we employ 
following iterative procedure [8] based on the Newton- 
Raphson method: 

Step 1 : Obtain a linear solution of Eq. (3.66) for a 
constant At = , At, and using only the linear 
part of flux boundary conditions with 
assumption u = U. 

With values of U at all nodes, find u, . 
material property values, and time step At 
at all nodes. Also compute derivatives of 
material properties required in evaluation of 
coefficients of the Jacobian matrix. 


Step 2 : 
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Step 3 : Update the Jacobian matrix and the residual 
vector, and find a new solution. 

Step 4 : Check if specified tolerance (s) is (are) 
achieved. If yes, then go to next step. 

Otherwise, go to step 2. 

Step 5 : With converged solution at current time-step, 
update initial conditions, and carry out 
steps (1) - (4) for solution at next step, 
if desired. 

★ 

In the above, the residual vector refers to f (x ) 

★ 

where x is an approximate solution for nonlinear equations 

f(x) = 0 (3.68) 

The Jacobian or tangent matrix (J) is defined as: 




■ax . 

j 


(3.69) 


chapter . 


We present a Newton-Raphson algorithm in next 


3 . 5 Energy Residual - A Physical 'Error- Estimate 

In this work, we do not make any attempt to present 
numerical error estimates which we consider to be outside 
its purview. In the following lines, we present a physical 
error estimate based on conservation of thermal energy. We 
check the input-output balance of thermal energy. The 
unbalance, herein referred to as the energy residual, can 
provide us with a fair estimate of accuracy and physical 
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acceptability of the computed approximate solution [25 J • 

The total thermal energy of the heat conducting 
system under consideration can be calculated at any 
time t m from: 


E m = /pc u m (x) dq (x) 

a 


(3.70) 


where p, c, and u m denote the mass density, the specific 
heat, and the temperature at point x at time t m . Net inflow 
of energy across the boundary r during the time period 
t m ~'*'< t < t m is given by 

t m 

E 1 ? = / dt (/ q dr) (3.71) 

m „ -i 

t 1 r 

q being the heat-influx, k(3u/3n). 

If rate of volumetric heat generation be q , then total 
energy generation in time-period t m t < t m is: 


E 


t 

= / 


m 


m-1 


dt 


(/ q d«) 

Q y 


(3.72) 


Assuming that all other forms of energy contained in the 
system remain constant during the interval t m 1 < t < t m , . 
conservation of energy requires: 

E m - (E m_1 + E m + E ) =0 (3.73) 

v m g 


where 
time t 


E m - 
m— 1 


1 


represents the thermal energy of the system at 
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Relation (3.73) is satisfied if solutions obtained were 
exact. With approximate solution the above relation 
is not satisfied exactly. We define the residual of 
the input-output balance by 

R m = E m - (E m_1 + E? + E^ ) (3.74) 

m g 

We call R 111 as the energy residual. We can clearly see 
that if R^ 1 = 0, then we have exact conservation of thermal 
energy. Furthermore, smaller the magnitude of R m , closer 
is the satisfaction of energy conservation, and this in 
turn, implies the accuracy of approximate solution. 

We may note that the above approach can be used in 
case of linear as well as nonlinear conduction problems. 

For evaluation of integrals in Eqs . (3.70) - (3.72) we 

employ finite-element approximations. 

We now proceed to the next chapter in which we discuss 
some important computational aspects related to numerical 
integration, solution of system of equations, approximations 
of variations of physical properties in case of nonlinear 
problems etc. We also present a macro-structure of the 
computer program developed. 



CHAPTER 4 


COMPUTATIONAL ASPECTS 


4 . 1 Evaluation of Boundary Integrals 

The integrals of following form need to be evaluated 
during the actual computations: 


I e (q, u ,£,) = f 0 S (x) u ( 5 , x) dr (x) (4.1) 

1 e 


I e (u, q*,£) = f T 0 S (x) q* (£,x) dr (x) (4.2) 


where Kernel functions u ( 5 , x) and q ( 5 , x) for two- 
dimensional problems are given by 


u*U,x) = . §7 In (7) 


* 

q U,x) 



(4.3) 


with d = ( 5 - x) .n , n being the unit outward normal, and 

r = | 5 -x| . 5 denotes the collocation point, x, the field 

point, and 0 e (x) represents the interpolation function. 

We note from Eq. (4.3) that Kernels u ( 5 , x) and q ( 5 , x) 
remain regular as long as the collocation point E, and the 
integration element r g do not coincide. Also that 
farther the integration element from the collocation point. 
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flatter is the variation of integration Kernels over it. 

'k k 

This nature of functions u and q suggests the use of 
selective numerical integration scheme. When the integra- 
tion element contains the collocation point or is very 
close to it, a standard numerical scheme fails to give 
accurate results. This situation warrants special attention. 

Approaches which have been suggested and used in 
the literature can be classified into four broad categories 
[33 ]. The first approach consists in the use of analytical 
integration over the elements which contain the collocation 
point or lie within a minimum distance from it, and, standard 
numerical quadrature over the rest of the elements. The 
analytical integration is, however, possible only for very 
simple elements. 

The second approach advocates the use of a special 
scheme - analytic ,semi~ analytic or numerical - for "the 
singular elements" (i.e. those elements to which the collo- 
cation point belongs) , and the use of a standard numerical 
scheme for the rest. This scheme has some problems when 
the collocation point is very close to the integration 
element . 

The third approach is to use special numerical 
integration schemes for the singular elements and also 
for the closest of the rest according to some specific 
criterion such as ratio of the distance between the 
collocation point and the integration element to the size 
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of this element, limitation of the truncation error of the 
integration etc. 


The fourth approach employs a standard numerical 
integration scheme with varying number of integration 
points over different elements (with increasing number of 
points of integration over the closest elements to the 
collocation point) according to one of the previous criteria. 
This is referred to as adaptive numerical integration. 

In present work, we use the analytic approach for 
evaluation of singular integrals, and a form of the fourth 
approach for the rest. In following lines we briefly 
outline the methods used. 


4.1.1 Regular Boundary Integrals 


The integrals over the elements not containing 
the collocation point are regular. We employ standard 
Gaussian quadrature with varying number of integration 
points such that farther the integration element lies 
from the collocation point, the fewer the number of 
integration points: 

1 n G 

f (y) dy = 

-1 k=l 


f fCy) dy = E 


(4.4) (a) 


where y k , w^. and n^ denote the point of integration, the 
weight corresponding to integration point y^., and total 
number of points of integration. We use the following 
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heuristic criterion suggested in references 10,34,35 based 

on the ratio (s) of distance of the integration element 

from the collocation point to length of the integration 

element for selection of number of integration points (n ) for 

G 

linear elements: 


s _< 1.5 

1.5 < s 1 2.5 

s > 5.5 


with 


s 


U- x c l 

1 



(4.4) (b) 


where x c is centroid of the integration element and 1 its length. 

For integration over elements very close to the collocation 
point such that the minimum distance of the integration element 
from the collocation point is less than the length of the element, 
we can divide the element into number of sub-elements, do the 
Integration over each sub-element using above method, and add 
the contributions of all sub-elements [ 33 ]. Division could be 
either uniform or graded. A criterion for uniform sub-division 
jiven in Reference [ 33] is : 


N . IF IX [l/d .1+1 (4.5) 

sub ' mm 

•/here d denotes the distance of collocation point from integration ele 
nent, IFIX denotes the integer part of the quotient, and 
the number of sub-divisions. A more sophisticated adaptive 
lumerical integration scheme and related references can be found in 
leference [33 ]. 
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4.1.2 Singular Boundary Integrals 

Over the element containing the collocation point, 

: k k 

the integral Kernels u and q become singular. The 

k 

integral containing u as its Kernel has integrable 

k 

(logarithmic) singularity, whereas the one having q as 

its Kernel can be evaluated only in Cauchy principal 

value sense. In what follows, we present a brief review 

and 

of methods available,/ results obtained with the one used 
in our present work^ 

A . Boundary Integrals with Integrable Singularity 
For evaluation of integral 

I (q, u* , 5) = / r 0 S (x) u* (£ , x) dT(x) (4.1) 

e e 

over singular element, we can use adaptive numerical 
integration, semi-analytic, or analytic integration, or 
special numerical integration schemes. Adaptive numerical 
integration schemes use a large number of integration 
points near the singularity. Semi -analytic methods 
make use of analytical integration over a short segment of 
integration element, and standard numerical quadrature 
over the rest. For singularity at hand, i.e. In (p) , 
we can make use of logarithmic weighting functions [36] 
in conjunction with Gaussian quadrature. For very _ simple 
elements, such as constant or linear elements, we can 
evaluate this integral analytically [ 10, 34, 35]. 
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To utilize Logarithmic Gaussian quadratures, we 
need to make some changes in integral (4.1), as these 
quadrature formulas apply to interval [ 0 , l] with singu- 
larity at zero. Following Pina [37], we put 


In ( 1/r) = In (x/r) + In (1/x) 


(4.6) 


Now, we can write: 

I e ^' u >0 =. I er (q, u*,£) + I es (q,u*, 5 ) ■ (4.7) 

with I = / - In (x /r) 0 e (x) J (x 1 ) dx 1 (4.8) 

er 0 2 tt e 

I = / ~—ln (1/x) 0 e (x) J (x 1 ) dx 1 (4.9) 

es z it e 


In the above J (x‘) denotes Jacobian of transformation 
e 

mapping r e onto interval [o,l] . Since x/r (x) remains 
bounded as x + the integral in (4.8) is regular, and 
standard Gaussian quadrature can be used for its evalua- 
ation. Singular integral in Eq. (4.9) has the proper 
form required by Logarithmic quadrature. 

The above approach is useful for any choice of 
boundary elements. However,, for constant and linear 
elements, we can easily use analytical integration 
procedure. For these elements, results of analytical 
integrations are : 
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A . 1 Constant Elements [lolfFia. 3.1(a)] 

I e (u*5 ) = - In (1/2) 1 (4.10) 

A . 2 Continuous Linear Elements [ 10 , 34 1 [Fig. 3.1(b) 3 

When collocation point is node 1 of the 
integration element: 


1 * 

Ig (u , 5) 


h [ 2~ ln(1):l 


(4.11) 



e 


* 

(u , £) 


m,i,] 


(4.12) 


Similarly, when collocation point coincides 
with local node 2 of the element, we have: 

Ig(u*,£) = jM^-lnU)] (4.13) 

I e (u *'5) = 47 -rf-lnd)] (4.14) 

A. 3 Partially Continuous Linear Elements [Fig. 3.1(c)] 
Elements of Type 1 : 

Case 1 : When collocation point coincides with 
freedom node 1 of the element: 

X e (u *'5 > = -d [ f - I lnU> ] C4 ‘ 151 

led, | ) = 47 - t f - f 1«U) 3 C4.161 



51 


CENTRAL L-RARY 

! i.T.. KA'»?'d* 


Acc. 


No. A. XO *4.2.1 A. 


Case 2 : Collocation point is local node 2 of 

the element: 


1 * 

T e (U } 


O * 

I e (u ,5) 


[<T 


— [(- 
4 it l v 3 


f ln3) 


ln3) - 


In (1/4) ] 


| In (1/4) ] 


(4.17) 


Elements of Type 2 ; 

Case 1 : Collocation point is local node 1 of 
the element: 


1 * 

I e lu ,5) 


9 * 

I e (U ,o 




4 / c ( l - f ln3 > 


In (1/4) ] 


| In (1/4) ] 


(4.18) 


Case 3 : Collocation point is local node 2 of 
the element: 


1 * 

Ig Cu 5) = 


9 * 

I 3 <u ,5) 


Xrl 
4 TT *“ 3 


In (1) ] 


4u ^3 3 ln (1) ^ 


(4.19) 


In all the above expressions, 1 denotes the length 


of the element. 
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B - Ca uchy Principal Value Integrals 
Integrals defined by ( 4 . 2 ), 

I e (u,q ,5) = / 0 e (x) q ( 5 , x) dr (x) (4.2) 

exist only in Cauchy principal value sense when collocation 

point lies on integration element. These integrals can be 

computed by 'using the finite parts of the integrals involved. 

Numerical integration formulas for finite part integrations 

are given in References [10,38] . For our present work, 

in which we use either constant or linear isoparametric 

* 

elements, over the singular element, q (£,x) is identically 
zero owing to orthogonality of vectors r and n. Thus, for 
constant and linear element: 

I e (u, q , ? ) = 0 (4.20) 


whenever the collocation point and integration element 
coincide. 

4 . 2 Evaluation of Domain Integrals 

Presence of heat generation (source/sink) inside 
the domain necessiates evaluation of domain integrals. 

When we use finite element type approximation for source 
term as discussed in Chapter 3, we are required to evaluate 
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integrals of the form: 

I k = ft' x k u * d ft (4.21) 

we again observe that this integral becomes singular when 
collocation point belongs to the domain element g , as u 
contains a ln(l/r) factor. In the following lines, we 
discuss the numerical methods used to evaluate these inte- 
grals. We restrict ourselves to linear triangular finite 
elements . 


4.2.1 Regular Domain Integrals 

If collocation point 5 does not belong to the 
integration domain element g , integral in (4.21) is 
regular. For triangular domain elements, we use quadrature 
scheme due to Hammer et al . given in Reference [ 1 0 ] : 


1 1 n . i i u 

/ ( f f(c 1 ,C , S 3 ) &z 1 )dr: 2 = "2 . Z f ^ 1 ’ ? 2 ' ^ W i 

> o i=l 

(4.22) 


where C-'s are triangular coordinates defined in 
1 

Fig. (3.3), and w.^ are weighting factors. 

Using (4.22), we can write: 


I k 


A 


e 


z f .C3I w i 

i=l 


(.4.2 3) 


with f = . x k u * 


*± 


We adopt the following criterion for the number of 
integration point on the basis of ratio of distance of 
collocation point from centroid of the element to the 
characteristic length of element: 


If s < 3, : n = 7 (Quintic scheme) 

(4.24) 

If s > 3 : n = 4 (Cubic scheme) , 

| x | 

with s = — , and l g = / 2 A q . x c 

e 

is the centroid of the element, 1 is its characteristic length, 

0 

and A is the area, 
e 


4.2.2 Singular Domain Integrals 

When collocation point C belongs to the domain 

element o , the integral 
0 

I* = / Y k u* drtx) (4.21)' 

becomes singular. To evaluate it, we follow the method 
proposed by Pina [ 37 ] . Let singular point coincides with 
the origin under isoparametric transformation. Letting 
r to denote the distance of the field point x from origin 
in isoparametric reference element, we can write: 

= f X V dfi = . / X k 2 ~“ln (1/r) J e (x)d (x) 

0 e U e 

= X v* J P 

g e x 2lT 

+ f x k ~2tT~ ln( - 1/: ^ J e (x) 


let. + le.s 


(4.25) 
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where Q denotes isoparametric reference element Fig. 4.1 . 
We choose a special reference element shown in 

A 

Fig. 4.1(c). Jacobian of transformation in 2 A . r 

e 

in (4.25) is defined by (with reference to Fig. 4.1 (c) ) 


/ A 2 
x 


+ y 


(4.26) 


We observe that first integral in right hand side of 
Eq. (4.25) is regular, and hence, method discussed in 
the previous section can be used for its evaluation. 

To perform integration using Hammer's quadrature scheme, 
we transform the integration domain into in 

natural triangular coordinates. Jacobian of transfor- 
mation is unity, and (x,y) and are related by: 


x = c , + S 2 / Y - 4 2 


(4.27) 


Therefore, 
I 


/ k __1 

* Y 


er a , 


x -^rln (r/r) J^ (x) dn (x) 


2 it' 


n 


A £ f ^2' ^'3^ ^ 2 . 

e i=l 


(4.28) 


where . 

f = 2 ^ — In (r/r), r = /x 2 +y 2 =’ / (c x + ^ 2 )2 + C 2 


; Isoparametric transformation 
: Special transformation 
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04 

X 


Fig.4.1 (a) Physical domain element. 

(b) Isoparametric reference element 

(c) Special reference element. 
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with 


X es = -r _[ (1+5 1> 27- ^a/t)x v a?! anj 

A e 1 2 w.w. f(?|, TV?) (4.32) 

/ 1 L 1 J 1 1 

i=l j=l 

f = . (1 +C 1 ) 2 ^— In (1/r) x k and 

r = . | (1 + 5 1 ) {1 + j (1 + r ^) 2 } 3 ^ 2 


Interpolation functions for linear approximation 
are defined in the terms of coordinates ( 5^, n^) by: 


X 


k 


x 1 = . \ (1 + 5 X ) (1 - V 

X 2 = 7 (1 + (1 + "i> 

X 3 = J (1 - 5 X ) 


( 4 . 33 ) 


In the procedure discussed above, we assumed that 
the collocation point is a vertex of the triangular 
element. However, even when the collocation point is 
on one of the sides or in the interior of the element, 
the domain element can be divided in two or three 
subelements, the integration procedure detailed above can 
be applied to evaluate integrals over each subelement and 
we finally add these to obtain the required integral. We 
employ this procedure of subdivision to evaluate integrals 
corresponding to the situation of collocation point belong- 
ing to one of the sides of the domain element, as will be 



the case if a constant or partially continuous linear 
element forms one of its sides. 
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A word about sequence of computation of integrals 
is in order. As we have noted that evaluation of an 
integral over an element requires computation of coordi- 
nates of integration points, its centroid etc., it is 
advisable that we take an element, and compute integrals 
over it from each collocation point. This will require 
lesser number of operations as compared to the sequence of 
choosing a collocation point, and computing integral over 
each element. We adopt same sequence of computation for 
boundary integrals also. 

4.3 Consideration of Symmetries 

By "symmetry" in this section, we mean planar 
symmetry with respect to global coordinate planes or 
sectorial symmetry. Symmetry conditions are used to 
reduce the number of elements needed to define the 
problem with consequent reduction in number of equations 
needed to solve it, and, hence, considerable reduction in 
computing time and memory requirements. 

Two approaches have been suggested and used in 
literature. The first one corresponds to the employment 
of knowledge of boundary conditions over the axes (planes) 
of symmetry, thereby, replacing the initial problem by 
a smaller one. This approach is quite general and can 
be applied to any kind of planar a sectorial symmetry. 
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The second approach consists in introducing symmetry 
conditions over the coordinates and the boundary conditions. 

We employ herein a direct condensation process with inte- 
gration over reflected elements. In this approach, no 
discretization of symmetry axes is required in direct contrast 
to the first approach. The number of equations required 
is, thus, reduced drastically. 

The second approach has the advantages of greater 
accuracy, and lesser number of equations. However, it 
requires greater number of integrations as compared to the 
former. Also, for each type of symmetry, a special consi- 
deration is required. If in addition, domain terms are also 
present, the process of discretization of domain and refle- 
ctions involved are not quite clear. In view of its res- 
tricted generality, we have adopted the first approach, 
which does not require any additional operations in program. 

4 . 4 Solution of System of Equations 

Depending on the nature of the diffusion process, 
and boundary conditions of the problem at hand, system of 
equations (3.48), (3.51) or (3.66) will be either linear 

or nonlinear. The system matrix is invariably full. Also 
that for 2-D problems, the size of system matrix is usually 
not very large as compared to the comparable system of 
equations obtained using Finite Element Method. Thug, in 
most of the cases, a direct solution procedure based on 
Gauss-elimination or its variants is usually the most suitable. 
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In our work, we have used a NAG Routine F04ATF based on Crout 1 s 
factorization method for solution of linear system of equations 

Solution of nonlinear system of equations requires a 
considerable effort in terms of selection of solution 
procedure, starting guess (as all of these methods are 
iterative in nature, and require a initial solution to start 
with) . In general, for weak nonlinear systems, modified 
Newton Methods are recommended, while for others, full 
Newton-Raphson Iteration is recommended. All these methods are 
strongly dependent on closeness of initial guess to the exact 
solution, and none of these has global convergence properly. 


4.4.1 N ewt on -R'a ph sou A'l qorithm 


Let, 


.f (x) = . 0 


(4.34) 


represents the nonlinear system of equations. The equality 

in above equation holds for exact solution vector x. To 

solve the above system we use iterative solution procedure. 

til 

Let x n be the approximate solution at n iteration. To 
obtain an expression for improved solution, let us expand f 

in Taylor series about x n . 


f (x) 


f tx n } 


+ c 


if ] 

3x J x=x. 


n 


(x-X ) + [ ] 

9x^ x=x 

n 

(x-x ) 2 + (4.35) 

n 
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Supposing that x is the exact solution, and 
neglecting terms containing second and higher order 
derivatives, we get: 


“ flx n> 


+ [-£] 


ax J x=x 


(x - 


n 


X n ) 


or. 


J A x 
n u n 


-f (x ) 
n 


(4 


with Ax = x - x , and J =[ l — ] is the Jacobian 

n n n L ax x=x 

n 

or tangent matrix. Having computed Ax, we can obtain 
the improved solution at (n+l) 1 " step by 


X n+1 


x + Ax 
n n 


(4 


We note that at each step, we have to solve a new set of 
linearized equations (4.36). 

The process of iteration is continued until the 
residual vector defined by f(x n+ ^) ^- s iciently 

small. 

To check the convergence of iterative procedure, 
we make use following vector norms for residual vector 
f and solution vector x: 

FNORM n = . Max [f^l , j = 1/...N ( 


.36) 


.37) 


.38) 
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XNORM 

n 


1 

N 


N 

I 

j = l 


n n-1 

x . - x . 

_J I 

, n . 

I x j I 


(4.39) 


where N represents the number of equations. Process of 
iteration is said to be converging if 


FNORM < 
n 


FNORM . 
n-1 


and 


XNORM < 
n 


XNORM 


n-1 


(4.40) 


If FTOL and XTOL denote the specified tolerances, then process 
of iteration ceases with message for converged solution if any 
or both of the following inequalities are satisfied: 


(i) FNORM < FTOL 

n (4.41) 

(ii) (FNORM < FNORM , and XNORM < XNORM n ) 
n n— 1 n n— 1 

and XNORM n < XTOL 

Process of iteration also terminates with an appropriate 
message if iteration is not making satisfactory progress. 

4.4.2 Initial Solution for Newton-Raphson Iteration 
The application of Newton-Raphson Algorithm 
discussed in above section requires a starting solution. 

In the present work, we compute a linear solution at each 
time-step of a linear system obtained by dropping the nonlinear 
part. We use this linear solution as initial solution for 
Newton-Raphson Iteration. 
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For problems with linear conduction and nonlinear 
boundary conditions, linear system is obtained by neglecting 
the nonlinear part of boundary conditions. For nonlinear 
heat diffusion problems, to obtain linear system for 
computation of initial solution we: 

(i) assume Atj = At for all j , 

(ii) assume u = U in flux boundary condition, and 

drop the nonlinear part of modified expression. 

The resulting system is linear in transformed variable U. 

4 . 5 Solution Procedure for Linear Conduction 

For physical problems with constant material proper- 
ties, system of boundary element equations - Eq. (3.48) 
or Eq. (3.51) remains linear as long as the boundary 
conditions are linear. We introduce prescribed boundary 
conditions in system of equations (3.48) or (3.51). 

4.5.1 Linear Systems 

For linear systems, we rearrange Eq. (3.48) or 
Eq. (3.51) after introducing the boundary conditions such 
that left hand side contains the unknowns, and right hand 
side is known. Resulting system is solved using a linear 
equation solver based on direct elimination or iterative 
procedures for fully populated real coefficient matrix 
system. We use NAG Library Routine FQ4ATF based on 


Crout's Method. 
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4.5.2 Nonlinear Systems 

Solution procedure for nonlinear systems involves 
two steps: 

Step 1 : Obtain a linear solution, which is used 
as the initial solution for step 2. 

Step 2 : Using linear solution from step 1 as ' 

initial guess for solution vector, use 
a routine based on Newton-Raphson Method 
detailed in Section 4.4.1. 

To obtain linearized system for step 1 above we 
use the procedure given in Section 4.4.2. 

4 • 6 Solution Procedure for Nonlinear Conduction 

Nonlinear heat conduction problems may be classified 
into two distinct categories: 

(i) Problems involving constant diffusivity, . 

(ii) Problems involving temperature dependent diffusivity. 

In the first category of problems, although the 
conductivity is temperature dependent, yet the variations 
of k and pc with temperature are such that (these differ 
at the most by a multiplicative constant) diffusivity 
is constant. Thus, time variable in transformed 


g 


equation 


3U 

3 T 


2 

V U + q 


(3.63) 
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is independent of space coordinates. The system of equations 
(3.66) is, therefore, same as that obtained for linear 
conduction problems, and procedure detailed in Section 4.5 
for linear problems is applied for solution of problems 
of this category. 

For problems with temperature-dependent diffusivity, . 
we outlined the iterative solution scheme in Section 3.4.2. 

In the following lines, we briefly discuss the computational 
details involved in iterative solution scheme. 

4.6.1 Iterative Solution Scheme 

Step 1 : To obtain linear solution which is used 
as the starting solution for Newton- 
Raphson Iteration, we put ATj = At in 
Eqn . (3.6 6)., Flux boundary ■ conditions 

in transformed space are linearized by 
putting u = U, and retaining only the 
linear part. By introducing this linear 
part in the appropriate system of 
equations, solve the resulting linear 

equation system, to get the linear 
solution. 

Step 2 : With linear solution obtained from step 1, 
employ Newton-Raphson Algorithm to obtain 
final solution. Expression for residual 


vector is 
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f(X m+1 ) = (C + ©H) U m+1 - © GQ m+1 - [C - ( 1-©) h] U m 

-(1-9) GQ m - B [(1-9) Q™ + © Q 171 ^ 1 3 

y y 

(4.42) 


where X denotes vector of unknowns. The coefficient of 
the Jacobian matrix are defined by 


j.m+1 

ij 


3 f. (X m+1 ) 

T xf 1 


(4.43) 


From (4.42), we have: 


, . « , v m+ 1 _m+ 1 

(l) when Xj = Qj 


J. . = - 9 G. . 

U ij 


/ ' • \ x. v m+l TT m+l 
(ii) when Xj = U 


■ 3 Q 


J. . = 
XJ 


At 


(U 


3 

m+1 

j 


C . . + 9 H. . — 9 G . . 
U XJ XJ 


m+1 

1 


3U 


m+1 


- U™ ) 


3U 


m+1 


’( Ax j ) 


j 

-1 


+ C . . 
XJ 


(4.44) 


with 


3 ( At .j ) 1 

3 U m+1 

J 




k 3 At 
J 


d k i 


d u m+1 
J 


C . 
J 


d C 


, m+1 
du . 

3 


d.p 

_L 


du 


m+1 


(4.45) 


8Q- 


m+1 


3U 


m+1 


k . 
J 


3 f c - (u, x, t ) 

3 U. 


and 
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4.6.2 Approximation of Property Variation s 

In the present work, we employ piecewise linear, 
approximation to model the temperature dependence of 
material properties. Derivatives of k, p , and c with 
respect to u, at any point j, are, thus, given by the 
slope of segment of curves to which Uj belongs. With 
reference to Fig. 4.2, the property value at temperature 
u lying on the segment i can be computed by: 


A(u) = A. + a. (u - u . ) 
11 1 


(4.46) 


where, a^ denotes the slope of segment 'i. 


To compute initial conditions and temperature 
boundary conditions in the Kirchhoff transform space with 
above variation for thermal conductivity, we first carry 
out a search to determine to which segment does a given u 
belong . Thereafter, we can write 


u 


U = T (u) = / k(u) du = 

u 


u . 

z 1 

U 


u 


k(x)d\ + / kU)dA 


u . 

l 


(4.47) 


(k, +k) 


U 


u. + 

i 


(u 


v 


U . 

with U = r 1 k(x) dx, assumed to be known. 

1 u 

r 

To perform an inverse transform, we again localize 
the segment i corresponding to- given value of U. Using 
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A 



Fig. 4.2 Piecewise -linear approximation of property 
(A) - temperature (u) curve. 
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Eq. (4.46) and Eq. (4.47), we obtain following expression 
for u 

1 -k. .+. ■ k 2 + 2a . (U - U. ) 

u = T 1 (U) = vu + - (4* 48) 

a . 
l 

j 

4.7 Program Structure 

In the following, we briefly discuss the macro- 
structure and organization of the Boundary Element 
program developed herein. The general scheme of the 
program is composed of a Main Program, and several 
Procedure Modules which can be run independently. The 
Main Program defines the dimensions of arrays and 
the common areas, and calls the rest of routines during 
the running of the program. Procedure Modules are 
usually composed of a general routine, which may control 
a set of other subroutines either specific for this 
module or shared with others. Each Procedure Module 
performs a specific function. Figure 4.3 shows the 
typical scheme for the program. 

In accordance with a scheme in Fig. 4.3, we have 
written a program which capable of solving a steady state 
or transient heat conduction problem in two-dimensions. 

The program can account for time-dependent boundary conditions 
and heat generation. It is composed of about 50 routines 
and 3000 lines. General flow chart of the program is shown 
in Fig. 4.4. Appendix contains FORTRAN listing of the 
Main Program "DRBEM" . 
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Fig. 4.3 Scheme of the programe. 




4.4 GENERAL FLOW CHART OF THE PROGRAM 


START 


Definition of arrays 
and COMMON areas 


nput .generation and 
treatment of data 


Set up integration points 
for Gaussian and Hammer’s 
quadratures 


Computation of 
domain integrals 


See if heat generation 
is present 


Computation of coordi- 
nate function matrix 
.related matrices 


Give time-increment 
and number of steps 


Compute b.c.’s in Kirc- 
hhoff transform space 
or heat generation term 

i * . A * a. ’ 


See if problem involves non 
linear conduction or time- 
dependent b.c.’s or source 
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i A ■ . 


,B. . 



FIG. 4.4 GENERAL FLOW CHART OF THE PROGRAM 




CHAPTER 5 


RESULTS AND DISCUSSIONS 


We present in the following sections results of analyses 
performed with the program DRBEM (outline of which has been 
presented in Chapter 4) . We have checked this program using 
three test-problems representing three possible cases of 
linear conduction, nonlinear conduction with constant diffusivity 
and nonlinear conduction with variable diffusivity. Section 5.1 
contains results of these test runs. Having checked the working 
of the program, we next present the results of studies on the 
convergence of the method with refining discretizations (keeping 
time-step fixed) , and with decreasing time-step (keeping spatial 
discretization fixed) in Section 5.2. Section 5.3 brings 
results of comparative performance of various time-stepping 
schemes. Finally, in Section 5.4, we present numerical solutions 
for nonlinear problems. 


5 . 1 Testing of the Program 
Test-Problem 1 : 

This non-dimensional problem, used in References C 40, 4 13 , 
is defined by 


3u 

5T 


3 2 u 

3X 2 


0 ■ in 0 < x^< 1 


subject to the initial and the boundary conditions: 


. (5.1) 
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u (x, 0) = sin ttx^ (5.2) 

u(0,t) = u ( 1 , t ) = 0 (5.3) 

Two-dimensional model of this problem together with appropriate 
initial and boundary conditions is shown in Fig. 5.1. The 
model has accounted for symmetry about x = 0.5. The exact 
solution is given by 

u (x ,t) = exp (— ir t ) sin x^ (5.4) 

Table 5.1 shows calculated boundary temperatures along x^-axis 
at two time instants obtained with time-step A t = 0.01, 
and, 20 equal linear boundary elements. Results are, in 
general, within 3% of the exact solutions and, thus, attest 
the correct working of the program for linear conduction 
problems. 

Test' Problem 2 : Nonlinear Bar 

This one-dimensional problem of a long bar with initial 
temperature of 0°C, subject to a unit heat input at the face 
x = 0 and insulated at x = L has been studied in References 
[8,42] . The geometrical dimensions and physical properties 
has been chosen as: L = 3m, K =(1.0 + 0.5 T)W m ^ °C 
c =(1.0 + 0.5T)J kg - ^ °C -1 and p= 1 kg m - ^. The bar was modelled 
as two-dimensional, with height 0.2m. This problem, representing 
the case of nonlinear diffusion with constant diffusivity 
was solved in Reference [8] using linear boundary elements 
of the size AT = 0.2m. We, instead, have modelled this 



° Geometric node 
x Freedom node 


* 


xc 

* g &- * * 

(b) 

Fig.5.1 Test problem-1 

(a) Two-dimensional model 

(b) Discretization 




Table 5.1. Boundary Temperatures along x.. -axis at Two Time Instants : Test Problem 
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problem using constant elements of the size AT = 0.3 m. 
Analysis has been performed in non-dimensional form, and, 
results are presented in Table 5.2 in dimensional variables. 
The exact solution at x^ = 0 is given by: 

T = 2'{[l + 2 (t / tt ) 1/2 ] 1/2 _ ± } (5 . 6) 

where T is the temperature in °c and t is the time in seconds. 

Results show a large error at first step which 
decreases to 0.1% at third step, and, then again shows an 
increasing trend. The initial error is due to the fact in 
numerical solutions ,unit flux is applied linearly over the 
first time-step. The latter trend can be attributed to the 
coarseness of the discretization and interpolation employed. 

We can still observe that errors are well within 15%, which 
is quite satisfactory in view of the discretization employed. 
The results of this problem further attest the sound 
working of the developed program. 

Test Problem 3 : Nonlinear Wall 

This problem has been taken from References [8,42 ] 
to check the performance of the program for general nonlinear 
conduction problems. The wall is defined by the domain, . 
n = { (x 1 , x 2 ) | 0 <_ x x £ 20, . 0 <. x 2 < 1 } . Initial 

and boundary conditions are: 



Table 5.2. Temperature Variation (°C) of face x = 0 of a Nonlinear Bar : Test Problem 

At = 0.05 s 


CM 

! 



vO 
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Results from Reference [ 8 J, Table 
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T 

(x 1 , x , 0) 

= 100 

(5.7) 

T 

(0, x 2 , t) 

= 200 for 0< t <10 




= 100 for t > 10 

(5.7) 

T 

(20, x 2 , t) 

= 100 


q 

(x 1 , 0, t) 

= q (x^, 1, t) = 0.0 



In the above x^ and are geometrical dimensions in cm, 
t is the time in seconds, and T is temperature in °c . 
q denotes heat flux in W cm . Material properties are 
conductivity K = 2+0.01T (W cm ^ °C and heat capacity 
pc = 9 (J cm -3 °C _1 ) . Taking T = 0°C as the reference 
temperature, and length of the wall as the characteristic 
length, we have non-dimensionalized the above variables, 
parameters and specified functions. Analysis has been 
performed in non-dimensional form. 

Results of analysis are presented in dimensional 
form at time t = 10s and time t = 13s, and, compared with 
finite element solution obtained with 20 equal linear 
finite elements [42], and boundary element solution using 
equal quadratic elements of length AT = 1 cm [ 8 ] in Table 5.3. 

Present analysis has used equal constant elements of size 
AT = 2 cm. We clearly observe that despite the crude 
discretization employed, results are fairly close to the 
reference solutions. This observation shows the correct 
working of the program for general nonlinear conduction 
problems also. 
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5 . 2 Numerical Convergence of Boundary Solutions 

Having tested the program DRBEM, we have next 
performed numerical experiments to study the numerical 
convergence of boundary solutions with refining discre- 
tization and decreasing time-increment. In following 
lines, we discuss the results obtained with fully- 
implicit (© = 1) and least-squares time integration 
schemes. The following linear problem has been 
considered: 


3 u 

at 



0 in Q = { x : | x^ | < 1 , | x 2 | < 1}, x= (x^x^ 

(5.10) 


with initial and boundary conditions: 


u (x, 0) = 1 on 0 (5.11) 

u (x,t) = 0 on r u x (5.12) 

Owing to the symmetry of the problem, we have solved 
Eq. (5.10) on a quadrant = (x j 0 <_ x^<_ 1, . 

0 <_ X 2 1 with appropriate initial and boundary 

conditions shown in Figure 5.2(a). The exact solution 
from Carslaw and Jaeger [43] is given by 
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u “ tt ^ ^n 1 -!- 1 GOS ttx ) exp [-(2n+l) 2 ] 

n=0 z i ft 


4 00 / 1 \ ^ o n j i i 

x tT ^ 2n + 1 cos tt x 2 ) exp [- (2n+l) 


2 it 


4 

(5.13) 

For small t, rapidly converging solution is given by [43] 


] 


» (2n+l)-x 

u = [ 1 - i (-1) { erfc - 

n=0 2/t 


+ erfc 


(2n+l) + x 1 
— 


} ] 


(5.14) 


x [1 - e (-l) n {erfc 
n=0 


(2n+l) -x 2 
2/t 


+ erfc 


(2n.+.l.) +x 0 

— 27t- > 5 


(5.15) 


where erfc (v) is the complimentary error function. 


5.2.1 Convergence with Mesh-Refinement 

For boundary element solution, we first employed 
a crude grid of 8 boundary elements with element size 
at = 0.5. To observe the convergence of boundary solutions 
with mesh refinement, we further subdivided the original 
mesh much more such as : at = 0.25 (16 element) and AT= 0.125 

(32 elements) for subsequent computations. Figs . 5 . 2 (b) - (d) 
show these meshes. 
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With f ully-implicit (© = 1) point-collocation 
scheme, we kept fixed time increment At = 1/128 for 
all computations. Boundary temperatures along x^-axis 
at the constant time t = 1.0 are plotted for the 
first two subdivisions in Fig. 5. 3. We can numerically 
observe uniform convergence of computed solutions to 
the exact solutions. Figure 5.4 shows the plot of 
relative error vs. boundary element size for the chosen 
subdivisions. Relative errors have been calculated 
at the point (0.1, 0.1). Slope of the regression line 
is 1:2, suggesting the second order of convergence with 
respect to mesh-refinement. 

* Jt , 

With the least- squares scheme, we have chosen 
At = 1/64. Figure 5.5 shows the boundary temperatures 
along x^ - axis for various subdivisions at time t = 1.0. 
We observe that the approximate solutions converge to 
the exact values, but the convergence is not uniform. 

To assess the convergence property of the 
approximate heat-fluxes, computed heat fluxes along the 
boundary x 1 = 1.0, 0 < x 2 < 1.0, at time t =1.0 are 

plotted for various discretizations in Fig. 5.6 for 
the fully-implicit scheme, and in Fig. 5. 7, for the 
least squares scheme. Time increments were At = 1/128 
for the fully-implicit scheme, and At = 1/64 for the 
least-squares scheme. From Fig. 5.6, we observe uniform 
convergence of boundary fluxes to the exact value 
with the fully implicit time integration scheme. With 






Temperature 



Fig. 5. 5 Boundary temperatures along x, - axis for various 
spatial grid sizes ( Least -squares scheme) 





0.6 0.8 


the boundary :x 
sizes( Point co 



Fig. 5.7 Heat fluxes along boundary x, = 1.0,0^ x 2 ^1 for 
various element sizes. (Least-squares scheme) 
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the least- squares scheme (Fig. 5.7), we once again see 
that the calculated fluxes converge to the exact values, 
but the convergence is not uniform. 

5*2.2 Convergence' with Decreasing Time -increment 

j 

To investigate the convergence of the approximate 
solution to the exact solution for the temperature and heat 
flux, we varied the time-increment while keeping the spatial 
mesh fixed. We chose the finest mesh shown in Fig. 5.2(d). 
Calculated temperatures along x^-axis are plotted in Fig. 5. 8 
and Fig. 5.9 for the fully implicit and the least-squares 
schemes respectively, for various time-increments. Fig. 5. 8 
shows uniform convergence of calculated boundary solutions to 
the exact values with the fully implicit time-integration 
scheme. From Fig. 5.9 we can observe that the solutions obta- 
ined using the least-squares scheme show oscillatory behaviour 
with a large time-step of At = 1/8. With the smaller time 
increments, however, we observe the uniform convergence of 
the boundary solutions to the exact ones with this scheme 
also. 

To assess the order of convergence with decreasing 
time-increment for the two schemes, we calculated relative 
errors at point (0.03, 0.0) for various time-increments. 

The slope of ..regression line in Fig. 5.10, for the fully 
implicit scheme, is 1:0.8, which suggests the first order 
of convergence. From the same figure, we observe that, 
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with the least- squares time-integration scheme, the 
regression line is a curve with non-constant slope. For 
range of time increments At = 1/16 to 1/32, the slope 
of regression line is approximately 1:2, suggesting 
the second order of convergence. Slope of the second 
segment of the regression line corresponding to the time- 
increment range of 1/32 to 1/64 is approximately 1:4. 

This suggests the fourth order of- convergence with 
respect to time-increments. These observations clearly 
indicate increasing order of convergence for the least- 
squares scheme with decreasing t ime- increment . 

Figures 5.11 - 5.12 show the plots of the computed 
heat fluxes along the boundary x^ = 0, 0 < ^ £ 1.0 for 
various time increments with the fully implicit and the 
least-squares schemes respectively. The results in 
Fig. 5.11 corresponding to the fully-implicit scheme 
indicate uniform convergence of the boundary fluxes to 
the exact solution. On the other hand, results obtained 
with the least-squares scheme (Fig. 5.12) are oscillatory 
for large time step, and converge to the exact solution 
with decreasing time step, but uniform convergence is 


not indicated. 



Temperature u x 10 



Fig. 5.8 Boundary temperatures along x, -axis for various 
time increments (Point collocation: 8=1 ) 


Temperature u x 10 



Fig 5 9 Boundary tempcraturesalong x, -axis for various 
time- increments (Least -squares scheme) 



Relative error 



Fig 5 10 Rate of convergence of the temperature (xi = 0.03,X2=i0) 



Fig. 5.11 Heat fluxes at the boundary :Xj =1.0, 0^x 2 ^1 for 

various time - increments (Weighted residual scheme:0=1) 




Heat f(ux-q x 10 



Fig. 5.12 Heat fluxes along the boundary x, = i 0 0<x <rin 
for various time increments ( LeaV-lquarfs ilhen 
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5 . 3 Comparison of the Least -Squares Scheme with other 
Weighted-Res idu'al' Schemes 

To compare the performance of the proposed least- 
squares scheme with the point collocation and the 
Galerkin schemes, the problem in Section 5.2 has been 
considered. The finest-grid in Fig. 5.2(d) has been 
employed. Calculated temperatures along x^-axis, at 
time t =1.0 with a time increment At= 1/64, are plotted 
in Fig. 5.13 for the fully-implicit (9 = 1), the Galerkin 
(9 = 2/3) , the Crank-Nicholson (9 = 1/2) , and the least- 
squares schemes. The least-squares scheme shows the 
greatest accuracy, but with associated larger-computational 
cost (computation time required with the least-squares 
scheme has been ^51 s as compared to ^46 s with other 
schemes for 64 time steps) . The Crank-Nicholson, the Galerkin 
and the fully-implicit schemes follow in order of decreasing 
solution accuracy 

To compare the computational cost associated with 
various schemes (which is directly proportional to 
the computation time), we have presented computation time requ- 
ired for various number of time-steps with the least- 
squares and the fully-implicit schemes in Table 5.4 
(Computation time with the Crank-Nicholson or the Galerkin 
scheme is almost same as that with the fully implicit scheme) . 
We can observe that the least-squares scheme requires larger 
computation time. However, the relative difference m 
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required time decreases with increasing number of time- 
steps. ( it is just 6% for 200 steps as compared to 35% 
for 10 steps) . Thus, in view of its higher accuracy, 
the least-squares scheme promises higher computational 
e fff c i enc y for large number of time steps. When number 
of time-steps is small, the difference in computation 
cost is marginal, and, a scheme offering higher accuracy is 
preferred. So, in this case too, least-squares scheme 
seems to be a better choice. 

Further, a comparison of Fig. 5.8 and Fig. 5.9 
reveals that the least-squares scheme obtains more accurate 
solutions with a time increment twice as large as that used 
with the fully implicit scheme. This observation further 
indicates the desirability of the least-squares scheme 
over the fully implicit scheme. 

Further, from Fig. 5.10, we observe that the rate of 
convergence for the fully-implicit scheme is of the first 
order. In comparison, the rate of convergence for the 
least-squares scheme is quadratic or higher. 

However, we note from Fig. 5.9 that least-squares 
scheme has produced oscillatory solution with large time- 
increment. Further investigation is called for under- 
standing this behaviour fully. At present, we add that 
with the least-squares scheme, caution must be exercised 
in selection of time— increment . 



Temperature u x 10 


100 



Fig 513 Boundary temperatures along x,-axis for various 

time- integration schemes. 




Table 5.4. Comparison of Computation Time for Various Number of Time-steps 
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5 . 4 Application to Problems with Nonlinear Boundary 
Conditions 


We consider an infinite slab of width 2L subjected 
to simultaneous convection and radiation at its surfaces. 
Dimensionless boundary condition has the form: 


with 


q = 


Bi 

Gr 


k — = - Bi (u + a )-Gr { (u+a) 4 - 

3 n c 


f (u) 
c 



' hL 


T 

a 

= if 

a T 

x b 

- T 

a 

e a (T - T ) 3 L 
b a __ 

T 

a 

- T 

c 

it 

o 

a c T 

b 

- T 
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4 


> 


(5.16) 


(5.17) 
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T-T 



a 


k 


K(T) 

K 

o 


where h is the convective heat transfer coefficient, K q 
the thermal conductivity at chosen temperature T q , e 
the relative emissivity, a. the Stef an-Boltzmann constant, 

T the temperature of convecting fluid, the temper- 

ature of the radiation source, and T a and T fe are suitably 
chosen reference temperatures for non-dimensionalization . 
Bi denotes the Biot number, and Gr is referred to as 
the radiation parameter. We have taken Bi = 4.0, and 
Gr * 4 _o, initial temperature of the slab has been 
taken to be T q = 1000. T c and T r have been assumed to 

be zero. We have chosen = 1000 and T a = 0 . Thus, 

= 0, and the initial condition 


we have a 


a_ * a. 
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u (x) =1.0. (5.18) 

o 

Owing to the symmetry with respect to the mid- 
plane of the slab, only half portion of width L needs 
to be analyzed. Figure 5.14(a) shows the geometric 
parameters, initial and boundary conditions, and two- 
dimensional rectangular region chosen for analysis. 

This two dimensional model has a length of 1.0 and 
height of 0.2. The boundary has been discretized using 
24 linear boundary elements. For computation of energy 
residual using finite element type approximation, the 
problem domain has been discretized using 40 linear 
triangular cells. 8 internal nodes in addition to the 
one corresponding to constant function have been chosen. 

Figure 5.14(b) shows the discretization employed. 

For both of following problems, the fully implicit 
time-integration scheme (9 = 1) with time— increment 
At = 1/128 has been employed. 

5.4.1 ' Linear Slab 

For linear slab, material properties are independent 
of temperature. In this case, dimensionalless properties are all 
equal to unity. Table 5.5 summarizes the results of the 
analysis in which temperature values at two extreme nodes 
node A (1.0, .0.1) and node B (0.0, .0.1) • - are presented. 

The calculated temperatures (denoted as DRBEM) are compared 
with the boundary element solutions using time-dependent 
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fundamental solution (denoted as BEM) by Onishi and 
Kuroki [25] in Fig. 5.15. The comparison shows quite 
good agreement between the two solutions. 

Furthermore, from last column of Table 5.5, we 
can see that the energy residuals are very small (less 
than 0.7%}. This clearly shows the satisfactory main- 
tenance of input-output energy balance, which, in turn, 
suggests that the computed solutions are accurate. 

In the above analysis, the Newton-Raphson algorithm 

detailed in Chapter 4 has been utilized. With a tolerance 
_ 4 

of 10 for relative norm of increments of unknown vector, 
number of Newton-Raphson iterations was 2. 

5.4.2 Nonlinear Slab 

For nonlinear slab, following variation of properties 
has been assumed: 

k (u) = 1.0 - 0.3 u 

p (u) = 1.0 (5.19) 

c (u) = . 1.0 

Results of the analysis are presented in Table 5.6. 

Figure 5.16 shows variation of temperatures at two extreme 
nodes - node A (1.0, 0.1) and node B (0.0, 0.1) with time. 

A tolerance of 10' 4 for the relative norm of the increments 
has been specified. Average number of Newton-Raphson 
iterations has been 4. 
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From the last column of Table 5 . 6, we can observe 
that relative energy residual is less than 1%. The input-' 
output energy balance is, thus, maintained satisfactorily. 
This, in turn, suggests high accuracy and physical accep- 
tability of the approximate solution obtained. 

From Fig. 5.16, we observe the same tendency as in 
Fig. 5.15 for linear slab with respect to the transient 
behaviour of boundary solutions. This observation further 
reinforces the physical acceptability of computed solu- 


tions . 


Simultaneous convection 



o Geometric nodes (Boundary) 
x Freedom nodes 
• Internal points 

Fig. 5. 14 (a) Geometry, initial and boundary 
conditions. 

(b) Discretization. 
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Table 5.5. Calculated Temperatures for Linear Slab subjected to Simultaneous Boundary 
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CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 

6 . 1 Conclusions 

The results of the present study have shown that J 
our attempt to develop a computer program based on the 
dual reciprocity boundary element method for general 
transient heat conduction problems in two dimensional 
isotropic medium has been successful. We have extended 
the method to deal with problems involving nonlinear 
boundary conditions. We have successfully incorporated 
the least-squares time integration scheme for linear 
problems. A systematic numerical study of convergence 
behaviour of the numerical solutions has been made. The 
following are the specific conclusions of the present 
work: 

( j_ ) The program developed herein has successfully 
implemented the dual reciprocity boundary 
element method to solve transient nonlinear heat 
conduction problems - involving nonlinear material 
as well nonlinear boundary conditions. The program 
provides the choice of constant and linear boundary 
elements for spatial discretization. It provides 
£qj~ the weighted— residual as well as the least 
squares time integration schemes (the latter for 
linear problems only) . The program has successfully 
accounted for time— varying boundary conditions. 




112 


Performance of the program with respect to the 
problems involving heat generation cannot be 
commented on, as no tests with such problems 
have been performed. The program, thus, provides 
many features of a comprehensive program. 

(ii) As mentioned above, a least-squares scheme has 
been successfully incorporated for linear problems. 
This has been done in the context of two point 
time integration schemes. 

(iii) From numerical experiments, with the f ully-implicit 
time integration scheme, the second order uniform 
convergence of boundary solutions with respect to 
mesh-refinement is obtained. With respect to 
time-increment, the above scheme shows first order 
convergence of the boundary solutions. 

With least-squares scheme, the convergence 
of the boundary solutions with refining spatial grid 
is not uniform. With respect to time-increment, 
this scheme shows atleast quadratic rate of conver- 
gence. With very large time-increment, this scheme 
has resulted in oscillatory solutions. 

(iv) As compared to the fully-implicit, the Galerkin 
and the Crank-Nieholson schemes, the least-squares 
scheme has been shown to be a computationally more 
efficient alternative in the context of linear 


problems . 
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(v) Problems involving nonlinear boundary conditions 
have been successfully solved. Accuracy of the 
computed solutions has been suggested by satis- 
factory maintenance of the input-output energy 
balance. 

6 . 2 Scope and Limitations 

The program developed herein can be used to analyze 
the thermal transients in the practical problems involving 2-D 
single homogeneous isotropic region. With nominal additions 
and modifications, this program can be used to analyze 
three-dimensional axisymmetric problems as well. 

The program, however, cannot be applied to multi- 
region problems in its present form. Also, the limitations 
imposed by fixed-storage assignation will make it difficult 
to run this program for large number of boundary nodes on 
a medium size machine. 

With its modular nature, the program offers immense 
flexibility in terms of modifications and enlargements. 


6 . 3 Recommendations for Further Study. 

The computer program developed herein requires 
specific attention with respect to the storage and 
management of the data. To handle problem— discretizations 
with large number of boundary nodes, the optimum use of 

and use of auxiliary storage area by means 


in-core storage 
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of variable number of buffer areas are called for. 
Incorporation of higher order interpolations in space 
and time are desirable. Theoretical and numerical 
investigations of convergence and stability of the 
method are needed. Extensive modifications are required 
for analyzing the coupled problems and the multi-region 
problems. Based on our experiences from the present work, 
we enlist in the following lines few suggestions towards 
the development of a comprehensive program based on the 
dual reciprocity boundary element method, and the develop- 
ment of the method itself. 

A. Programming Aspects Storage and Management of Data 

We suggest incorporation of the following 
data structure [33] to take into account the restri- 
ctions on the in-core storage, easy checking and 
access to each of the arrays involved and use of 
auxiliary memory: 

Dynamic assignation of the in-core storage 
area during the running of the program. 

Use of auxiliary memory by means of a variable 
number of buffer areas. 

Use of a set of routines for initialization of 
the global control area, creation of a specific 
array, its location in the working area and 
location of a particular element of an array and 
its replacement, if necessary. 
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B . Numerical Integration 

Incorporation of adaptive numerical integration 
schemes is suggested. In the present work, a heuristic 
criterion for selection of the number of integration 
points has been employed. Instead, use of a criterion 
based on limitation of the truncation error of integration 
or ratio of the distance of the collocation point from 
the integration element to the length of the element 
is recommended. For quadratic and higher order elements 
incorporation of a finite-part integration scheme for 
evaluation of Cauchy principal value integrals is 
required . 

C. Interpolation and Shape Functions 

To better model the geometry of the problem 
domain and variation of field variables, incorporation 
of quadratic and cubic elements, and interpolation and 
shape functions is suggested. 

D . Axi symmetric Problems 

The present 2-D program can be easily modified 
to account for three-dimensional axisymmetric problems as 
well. To affect this extension, incorporation of 
routines to compute integral coefficients and coordinate 
functions for axisymmetric problems is recommended. 
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E. T i me- Integration Sc h erne s 

In the context of two-point time integration 
schemes, the least-squares scheme has shown consi- 
derable promise for linear problems. Further exten- 
sions for linear conduction with nonlinear boundary 
conditions and nonlinear conduction problems are 
recommended. We recommend theoretical investigations 
of stability and convergence of this scheme. Compu- 
tational viability of the scheme for nonlinear problems 
also needs to be investigated. 

Use of three-point time integration schemes 
comes as the next choice. Also recommended is the 
theoretical as well as numerical experimental 
stability analysis of these schemes. To account for 
the abrupt discontinuity at the start of time-inte- 
grations introduced by rapid change of boundary 
conditions, incorporation of a procedure for smoothing 
the discontinuity is recommended. Such a procedure 
in the context of the finite element methods is detailed 
in Reference [32] , and we suggest for adoption of a 
similar procedure in our context. 

F . Theoretical and Experimental Analysis of Convergence 
of Boundary Solutions 

Theoretical analysis of approximation errors, 
stability and convergence of boundary solutions in the 
context of the dual reciprocity boundary element 
methods remains an open area. Estimation of approxi- 
mation errors forms the first pre-requisite for the 
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development of the adaptive grid techniques based 
on this method. 

Numerical experiments are recommended to 
investigate the convergence behaviour of the solu- 
tions with spatial and time grid refinements for 
each choice of spatial and time-interpolation. 

G. Adaptive Grid Techniques 

As a next step in the development of the 
method, attempts towards the development of adaptive 
schemes for the time-dependent problems are called 
for . 

H. Extensions to Multi-region and Anisotropic 
Problems 

For analysis of temperature field in composite 
bodies using the dual reciprocity boundary element 
method, the following step by step procedure is 
recommended: 

Step-1: Set up equations of the form of Eq. (3.48) 

or Eq. (3.66) for each region assuming it 
to be a separate one. To affect this 
process, compute the integral coefficient 
matrices and coordinate function matrices 
for each region. 
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Step-2: Apply the compatibility of physical temper- 
atures and equilibrium of heat fluxes at 
each interface. Imperfect contact can be 
modelled as a convective resistance. 

Step-3: Assemble the coupled system of equation. 

Step-4: Apply the boundary conditions, and follow 

the same set of solution procedure as for 
the single region problems. 

The above solution procedure will require a 
careful ordering of the coupled system of equations. 

A systematic assemblage process is called for at 

the Step-3 above. The solution of system of equations 

requires further considerations. 

Extensions to anisotropic problems should also 
be considered. 

I . Analysis of Coupled Problems 

For analysis of the practical problems involving 
coupling with temperature field in a conducting region, 
the program developed herein can be suitably augmented 
with the one(s) dealing with the other field variable (s) 
For problems involving weak coupling such as thermo- 
mechanical problems in analysis of thermal stresses, the 
results of analysis obtained with our program can be 
used as the temperature history for computation of 


stresses . 
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